Cosparse seismic data interpolation

Many modern seismic data interpolation and redatuming algorithms rely on the promotion of transform-domain sparsity for high-quality results. Amongst the large diversity of methods and different ways of realizing sparse reconstruction lies a central question that often goes unaddressed: is it better for the transform-domain sparsity to be achieved through explicit construction of sparse representations (e.g., by thresholding of small transform-domain coefficients), or by demanding that the algorithm return physical signals which produces sparse coefficients when hit with the forward transform? Recent results show that the two approaches give rise to different solutions when the transform is redundant, and that the latter approach imposes a whole new class of constraints related to where the forward transform produces zero coefficients. From this framework, a new reconstruction algorithm is proposed which may allow better reconstruction from subsampled signaled than what the sparsity assumption alone would predict. In this work we apply the new framework and algorithm to the case of seismic data interpolation under the curvelet domain, and show that it admits better reconstruction than some existing L1 sparsity-based methods derived from compressive sensing for a range of subsampling factors. Introduction The recognition of the role of sparsity in the representation of seismic data under redundant transform domains, such as windowed Fourier, contourlets, curvelets, etc., have led to a well-developed field of inversion-based non-parametric seismic data interpolation methods (Abma and Kabir, 2006; Herrmann and Hennenfent, 2008; Li et al., 2012), where missing seismic traces in a record can be recovered without a-priori knowledge of the subsurface. While much work has been done on the specifics of the transform and the algorithms used to promote sparse representations, an interesting question remains unaddressed: is it better for the transform-domain sparsity to be achieved through explicit construction of sparse representations (e.g., by thresholding of small transform-domain coefficients), or by demanding that the algorithm return physical signals which produces sparse coefficients when hit with the forward transform? This is related to the well-known “synthesis or analysis” problem in signal processing. Most works on this topic indadvertedly prefer one over the other without much discussion of their relative merits. Some recent studies have implied that the latter “analysis-based” approach poses a stronger condition on the solution through, and potentially achieves better recovery from signal subsampling than approaches based on explicit construction of sparse coefficients. The resulting theory of cosparsity provides new recovery guarantees and algorithms for the analysis problem. We apply this finding to the problem of curvelet-based seismic data interpolation, and demonstrate an improved interpolation result using an algorithm based on the cosparsity framework. Theory The goal of seismic interpolation is to invert the under-determined system, y = Rf+n where f is the fully-sampled seismic signal, y the observed data, n some zero-meaned noise and R is a restriction operator that removes unobserved traces. For the sake of simplicity, we assume in this work that the observations lie exactly on the reconstruction grid. To obtain a unique solution, additional regularization has to be imposed in a way that favors the recovery of the correct fully-sampled seismic record. Assuming that the full seismic signal f permits a sparse (or compressible) representation through a transform domain defined by a synthesis operator S (a linear mapping of coefficients to physical signal), compressive sensing theory shows that it is possible to reconstruct f from y without any additional information such as subsurface velocities, depending on the particular properties of R and S (Herrmann, 2010). This led to an approach called Curvelet Recovery by Sparsity-promoting Inversion (CRSI, Hennenfent, 2008), where S := CT is the inverse Curvelet transform, and the solution is obtained by solving f̃ = S · argmin x ‖x‖1 subject to ‖y−RSx‖2 ≤ σ , (1) where ‖x‖1 := ∑i |xi| is the `1-norm of the coefficients, and σ is some predetermined noise level. This optimization problem is often referred to as “`1 synthesis” signal reconstruction, as it relies the possibility of constructing a sparse model of the signal through synthesis operator. Alternatively, an “`1 analysis” problem may be formulated that relies on the analysis operator Ω of a transform domain (mapping physical signals to coefficients). Its analogue of CRSI defines Ω := C as the forward Curvelet transform, and solves f̃ = argmin f ‖Ωf‖1 subject to ‖y−Rf‖2 ≤ σ . (2) This formulation appeared recently in Li et al. (2012). Surprisingly, these two forms often appear as interchangeable approaches in literature, even though it is known that they lead to different solutions for most choices of Ω and S (Elad et al., 2007). In fact, we can only assume that (1) and (2) give the same solution when we impose an additional constraint x = ΩSx on (1). This is trivial when Ω and S are respectively the forward and inverse of an orthonormal transform (where the two approaches will equate by definition), but is not generally true when using redundant transforms. For curvelets, this would necessitate x = CCT x, which is not true for most choices of x. It is thus evident that the analysis problem poses more constraints on the solution compared to the synthesis problem, and the findings of Li et al. (2012) seem to indicate that there is a small uplift from using the analysis approach. 75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013 London, UK, 10-13 June 2013 The cosparsity model for analysis problems: Rigorous analysis of `1 synthesis, such as recovery guarantees, had been better-developed than `1 analysis due to its association with the standard form of compressive sensing problems. A significant development for analysis problems emerged recently form the insights of Nam et al. (2013), which investigated the signal recovery guarantees of f̃ = argmin f ‖Ωf‖0 subject to ‖y−Rf‖2 ≤ σ , (3) where the `0 term measures the absolute number of non-zero coefficients in Ωf. They found that by only looking at the zeros of Ωf, they can derive signal recovery guarantees from undersampled signals from the solution of (3), in a similar fashion to the role of sparsity in synthesis problems and compressive sensing. The authors coined the name of “cosparsity” of signal f for the number of zeros present in Ωf. In general, cosparsity is a more powerful condition than sparsity, as the zero-elements of Ωf specify a subspace from which f must be orthogonal to. Although it is well-known from compressive sensing that the `1-norm is a good heuristic for sparsity in signal recovery problems, it is not clear whether the same property holds true for cosparsity. The authors therefore proposed a greedy algorithm that attempts to solve equation 3, otherwise a combinatorial problem, called Greedy Analysis Pursuit (GAP). The goal of GAP is to sequentially remove one or more rows of Ω that form large inner-products with f, until the reduced operator ΩΛ, where Λ is the subset of row indices remaining from the full N number of rows, produces coefficients close to zero in a regularized least-squares problem minx ‖ΩΛx‖2 s.t. ‖y−Rx‖2 ≤ σ . The rate of removal of rows from Ω is controlled by a selection factor t. Notably the authors have reported that using GAP a large improvement in signal recovery from undersampling can be seen compared to solving the `1 analysis problem (using standard convex optimization solvers). In the next section, we apply the GAP approach to solving CRSI-type problems for seismic data interpolation. Algorithm 1: Greedy Analysis Pursuit (GAP) Result: Approximate solution f̃ to the `0 analysis problem (Equation 3) choose noise level σ , selection factor 0 < t ≤ 1, stopping coefficient size ρ initialize iteration counter k← 0, row index set for analysis operator Λ0← {1,2,3, . . . ,N} f̃0← argminx ‖Ωx‖2 subject to ‖y−Rx‖2 ≤ σ repeat obtain coefficients: x ← Ωf̃k find the indices of the largest entries of x: Γk+1 ← {i : |xi| ≥ t ·max j|x j|} cull the selected indices from the set of remaining rows: Λk+1 ← Λk\Γk+1 update solution: f̃k+1← argminx ‖ΩΛk+1x‖2 subject to ‖y−Rx‖2 ≤ σ k← k+1 until remaining coefficients close enough to zero: max j|x j| ≤ ρ Numerical examples We use synthetic experiments to demonstrate uplifts in curvelet-based interpolation accuracy using by GAP in comparison to `1-norm based methods. A single fully-sampled synthetic shot record with 256 traces spaced 15m apart is used as the ground truth model for f. Our observed traces are determined from a random jittered-sampling scheme (Hennenfent and Herrmann, 2008) out of these traces, with the restriction mask R constructed accordingly. Interpolate in the common-shot domain is more challenging than in the more typical common-offset domain, as the apex of the hyperbolic reflection events introduces significant curvature that is often regarded as difficult to interpolate without knowledge of velocity. We chose this to better highlight any differences between the reconstruction methods. We compared reconstruction results based on three different algorithms: the approximated solution to equation (3) using GAP, the `1 synthesis solution to equation (1) using SPGL1 (van den Berg and Friedlander, 2008), and the `1 analysis solution to equation (2) using the analysis mode of NESTA (Becker 75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013 London, UK, 10-13 June 2013 et al., 2011). For GAP, the stopping coefficient size ρ is set to 10−5 times the largest coefficient magnitude of ΩRT y, while the other algorithms are set to terminate at default conditions. Figure 1 shows the results obtained from these different methods from 50% of total traces observed. At areas of low curvature in the far offsets, all methods show successful reconstruction. However, the GAP reconstruction results clearly show much better reconstruction of of the apex compared to the other two methods. This uplift is corroborated by numerical calculations of the signal-to-noise ratio in Figure 2 where GAP clearly shows superior reconstruction over the `1-based methods up to 50% missing traces. It is interesting to note that although the NESTA result is also based on the analysis form of the reconstruction problem, its solution are more or less comparable to the `1 synthesis result. This demonstrates the value of choosing an algorithm that is designed for a cosparsity-based model, and suggests that the `1-norm may not be as good a cosparsity heuristic as it was for sparsity. The main drawback of GAP is the higher iteration count needed compared to the other methods, perhaps partly due to its simplicity. Here we find that the selection factor t has a large impact on the efficacy of the algorithm in terms of speed and solution quality. Lowering t eliminates more coefficients per iteration and accelerates convergence up to a point where the solution quality begins to suffer. Figure 3 summarizes the effect of t on the convergence of GAP for this example, as well as the comparative performance between the different algorithms. Although the convergence of GAP is slower than its counterparts, it is able to achieve a lower error in the end.